── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Adjuntando el paquete: 'dagitty'
The following object is masked from 'package:rio':
convert
Adjuntando el paquete: 'ggdag'
The following object is masked from 'package:stats':
filter
Adjuntando el paquete: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Cargando paquete requerido: Matrix
Adjuntando el paquete: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Adjuntando el paquete: 'lme4'
The following object is masked from 'package:rio':
factorize
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
Cargando paquete requerido: carData
Adjuntando el paquete: 'car'
The following object is masked from 'package:purrr':
some
The following object is masked from 'package:dplyr':
recode
Cargando paquete requerido: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Adjuntando el paquete: 'brms'
The following object is masked from 'package:glmmTMB':
lognormal
The following object is masked from 'package:lme4':
ngrps
The following object is masked from 'package:stats':
ar
This is bayesplot version 1.12.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Adjuntando el paquete: 'bayesplot'
The following object is masked from 'package:brms':
rhat
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
##
## Synth Package: Implements Synthetic Control Methods.
## See https://web.stanford.edu/~jhain/synthpage.html for additional information.
Warning: package 'weathercan' was built under R version 4.4.3
As of v0.7.2, the `normals` column in `stations()` reflects whether or not there
are *any* normals available (not just the most recent).
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
Registered S3 method overwritten by 'lfe':
method from
nobs.felm broom
Adjuntando el paquete: 'plm'
The following objects are masked from 'package:dplyr':
between, lag, lead
Adjuntando el paquete: 'PanelMatch'
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:stats':
weights
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Cargando paquete requerido: bsts
Cargando paquete requerido: BoomSpikeSlab
Cargando paquete requerido: Boom
Adjuntando el paquete: 'Boom'
The following objects are masked from 'package:brms':
ddirichlet, rdirichlet
The following object is masked from 'package:stats':
rWishart
Adjuntando el paquete: 'BoomSpikeSlab'
The following object is masked from 'package:stats':
knots
Cargando paquete requerido: zoo
Adjuntando el paquete: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Cargando paquete requerido: xts
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Adjuntando el paquete: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Adjuntando el paquete: 'bsts'
The following object is masked from 'package:BoomSpikeSlab':
SuggestBurn
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
Adjuntando el paquete: 'forecast'
The following object is masked from 'package:brms':
ma
Code
#special repository indicated or the packageif(!require(weathercan)){install.packages("weathercan", repos =c("https://ropensci.r-universe.dev", "https://cloud.r-project.org")); library(weathercan) }#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_if(!require(bpmn)){devtools::install_github("bergant/bpmn")}
Warning in readLines(paste0(gsub("f1/", "", getwd()), "/key.txt")): incomplete
final line found on 'H:/Mi unidad/PERSONAL
ANDRES/UCH_salud_publica/pasantia/f1/key.txt'
API key has been set for the current session.
Forecast with Time GPT
Code
#https://cran.r-project.org/web/packages/nixtlar/vignettes/prediction-intervals.html#https://cran.r-project.org/web/packages/nixtlar/vignettes/exogenous-variables.html#https://cran.r-project.org/web/packages/nixtlar/vignettes/cross-validation.html#https://cran.r-project.org/web/packages/nixtlar/vignettes/anomaly-detection.htmlcollisions_weather_corr_rect_synth_tr$Y_lag10 <-collisions_weather_corr_rect_synth_tr|>mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>ungroup()|>pull(Y_lag10) collisions_weather_corr_rect_synth_tr$Y_lag10 <-collisions_weather_corr_rect_synth_tr|>mutate(third_value =nth(na.omit(Y_lag10), 1),Y_lag10 =ifelse(is.na(Y_lag10), third_value, Y_lag10) )|> dplyr::select(Y_lag10)|>ungroup()|>pull(Y_lag10)df35 <- collisions_weather_corr_rect_synth_tr %>%arrange(tr_contr_sens, year.x, yday_corr) %>%# 1,2,3… days relative to racefilter(yday_corr<36) |>transmute(unique_id ="treated", # exactly this nameds =as.POSIXct(as.Date("2019-05-19") + (yday_corr -1)),y = rate_veh,mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,mean_min_temp_mean_lin= mean_min_temp_mean_lin,mean_max_temp_mean_lin= mean_max_temp_mean_lin,Y_lag10= Y_lag10 )df35_post <- collisions_weather_corr_rect_synth_tr %>%arrange(tr_contr_sens, year.x, yday_corr) %>%# 1,2,3… days relative to racefilter(yday_corr>=36, yday_corr<=max_time) |># take the first 30transmute(unique_id ="treated", # exactly this nameds =as.POSIXct(as.Date("2019-05-19") + (yday_corr -1)),y = rate_veh,mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,mean_min_temp_mean_lin= mean_min_temp_mean_lin,mean_max_temp_mean_lin= mean_max_temp_mean_lin,Y_lag10= Y_lag10 )df33 <- collisions_weather_corr_rect_synth_tr %>%arrange(tr_contr_sens, year.x, yday_corr) %>%# 1,2,3… days relative to racefilter(yday_corr<34) |>transmute(unique_id ="treated", # exactly this nameds =as.POSIXct(as.Date("2019-05-19") + (yday_corr -1)),y = rate_veh,mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,mean_min_temp_mean_lin= mean_min_temp_mean_lin,mean_max_temp_mean_lin= mean_max_temp_mean_lin,Y_lag10= Y_lag10 )df33_post <- collisions_weather_corr_rect_synth_tr %>%arrange(tr_contr_sens, year.x, yday_corr) %>%# 1,2,3… days relative to racefilter(yday_corr>=34, yday_corr<=max_time) |># take the first 30transmute(unique_id ="treated", # exactly this nameds =as.POSIXct(as.Date("2019-05-19") + (yday_corr -1)),y = rate_veh,mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin= mean_median_total_precip_median_lin,mean_min_temp_mean_lin= mean_min_temp_mean_lin,mean_max_temp_mean_lin= mean_max_temp_mean_lin,Y_lag10= Y_lag10 )fc <-nixtla_client_forecast( df35,h =7,level =c(95),freq ="D", # <— this fixes the infer_frequency errorfinetune_steps =10, #Number of steps used to finetune 'TimeGPT' in the new data.finetune_loss ="rmse")
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
fc33 <-nixtla_client_forecast(rbind(mutate(df33[1:2,], ds=ds-3), df33),h =7,level =c(95),freq ="D", # <— this fixes the infer_frequency errorfinetune_steps =10, #Number of steps used to finetune 'TimeGPT' in the new data.finetune_loss ="rmse")
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
p <-nixtla_client_plot(rbind(df35, df35_post), fc,max_insample_length =100)p33 <-nixtla_client_plot(rbind(mutate(df33[1:2,], ds=ds-3), df33, df33_post), fc33,max_insample_length =100)# 1) Banda de CI (layer 1)p$layers[[1]]$geom_params$fill <-"#FFC10740"# amarillo semitransparentep$layers[[1]]$geom_params$alpha <-0.4# 2) Línea histórica y pronóstico (layer 2)# p$layers[[2]]$aes_params$colour <- "steelblue"p$layers[[2]]$aes_params$linewidth <-1.0p33$layers[[1]]$aes_params$fill <-"#eb8f8f"# verde semitransparentep33$layers[[1]]$aes_params$alpha <-0.5# 1) Banda de CI (layer 1)p33$layers[[1]]$geom_params$fill <-"#FFC10740"# amarillo semitransparentep33$layers[[1]]$geom_params$alpha <-0.4# 2) Línea histórica y pronóstico (layer 2)# p$layers[[2]]$aes_params$colour <- "steelblue"p33$layers[[2]]$aes_params$linewidth <-1.0p33$layers[[1]]$aes_params$fill <-"#eb8f8f"# verde semitransparentep33$layers[[1]]$aes_params$alpha <-0.5t0 <-min(p$data$ds)t0_33 <-min(p33$data$ds)t0_posix <-min(p$data$ds)t0_33_posix <-min(p33$data$ds)xint35 <- t0_posix +days(34)t0 <-as.Date(min(p$data$ds))t0_33 <-as.Date(min(p33$data$ds))# 3) Aplica scale_x_datetime() para transformar POSIXct → numéricop +scale_color_manual(values =c(# y = "#FFD600", # amarillo intensoTimeGPT ="red"# violeta (DarkOrchid / BlueViolet) ), labels=c( TimeGPT ="Prediction (95% IC in gray area)" ) ) +scale_x_datetime(name ="Days of follow-up (35 days before race,\n7 days after the race)",limits =c(min(p$data$ds), max(p$data$ds)), breaks =seq(min(p$data$ds), max(p$data$ds), by ="7 days"),labels =function(x) {# convierto POSIXct a Date y resto t0as.integer(as.Date(x) - t0 +1) } ) +theme_classic(base_size =13) +theme(strip.text =element_blank(), # quita el label de facetastrip.background =element_blank(),legend.position="bottom" )+geom_vline(xintercept = xint35,linetype ="dashed",color ="black",linewidth =1 ) +labs(y="Highspeed collisions per 1M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates
Code
p33+scale_color_manual(values =c(# y = "#FFD600", # amarillo intensoTimeGPT ="red"# violeta (DarkOrchid / BlueViolet) ), labels=c( TimeGPT ="Prediction (95% IC in gray area)" ) ) +scale_x_datetime(name ="Days of follow-up (33 days before race,\n7 days after the race)",limits =c(min(p33$data$ds), max(p33$data$ds)), breaks =seq(min(p33$data$ds), max(p33$data$ds), by ="7 days"),labels =function(x) {# convierto POSIXct a Date y resto t0as.integer(as.Date(x) - t0_33 +1) } ) +theme_classic(base_size =13) +theme(strip.text =element_blank(), # quita el label de facetastrip.background =element_blank(),legend.position="bottom" )+geom_vline(xintercept = xint35-days(2),linetype ="dashed",color ="black",linewidth =1 ) +labs(y="Highspeed collisions per 1M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediciton on the treated (adjusting on min max temperature and lagged precipitiation (2 days) and lagged rates
Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, Y_lag10
Code
p33_qc <-nixtla_client_plot(rbind(mutate(df33_qc[1:2,], ds = ds -3), df33_qc, df33_post_qc), fc33_qc,max_insample_length =100)## ── Estética idéntica a tu original ───────────────────────────────────────────p33_qc$layers[[1]]$geom_params$fill <-"#FFC10740"p33_qc$layers[[1]]$geom_params$alpha <-0.4p33_qc$layers[[2]]$aes_params$linewidth <-1.0t0_33_qc <-as.Date(min(p33_qc$data$ds))xint33_qc <-min(p33_qc$data$ds) + lubridate::days(32)p33_qc +scale_color_manual(values =c(TimeGPT ="red"),labels =c(TimeGPT ="Prediction (95% IC in gray area)") ) +scale_x_datetime(name ="Days of follow-up (33 before, 7 after)",limits =c(min(p33_qc$data$ds), max(p33_qc$data$ds)),breaks =seq(min(p33_qc$data$ds), max(p33_qc$data$ds), by ="7 days"),labels =function(x) as.integer(as.Date(x) - t0_33_qc +1) ) +theme_classic(base_size =13) +theme(strip.text =element_blank(),strip.background =element_blank(),legend.position ="bottom" ) +geom_vline(xintercept = xint33_qc,linetype ="dashed",color ="black",linewidth =1 ) +labs(y ="High-speed collisions per 1 M vehicles")
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Prediction on the treated (-2 days before race, 7 days after — Québec)
Causal impact
We used the CausalImpact package, useful for time series data, inferring a counterfactual post race period based on pre-intervention components1.
Code
# Model 2ss2d <-list()# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TRENDss2d <-AddLocalLevel(ss2d,c(as.numeric(df35$y), rep(NA,7)) ) ## Add weekly seasonalss2d <-AddSeasonal(ss2d,c(as.numeric(df35$y), rep(NA,7)), nseasons=7, season.duration =1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑOss2d <-AddAutoAr(ss2d, y =c(as.numeric(df35$y), rep(NA,7)), lags =10) #NO PUEDO AREGAR AR1 CON POISSON# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).# Add external regressors explicitly (e.g., temperatures, precipitation)covariates <-cbind(# y = as.numeric(df35$y),mean_min_temp_mean_lin = df35$mean_min_temp_mean_lin,mean_max_temp_mean_lin = df35$mean_max_temp_mean_lin,mean_median_lag_2_prec_median_imp = df35$mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin = df35$mean_median_total_precip_median_lin,Y_lag10 = df35$Y_lag10)# Ensure names explicitlycolnames(covariates) <-c(#"y","mean_min_temp_mean_lin","mean_max_temp_mean_lin","mean_median_lag_2_prec_median_imp","mean_median_total_precip_median_lin","Y_lag10")covariates<-as.data.frame(covariates)covariates_future <-as.matrix(df35_post[, c("mean_min_temp_mean_lin","mean_max_temp_mean_lin","mean_median_lag_2_prec_median_imp","mean_median_total_precip_median_lin","Y_lag10")])covariates_full <-rbind(covariates, covariates_future)# Include regression termsss2d <-AddDynamicRegression(ss2d, covariates_full)model2d_ratio <-bsts(c(as.numeric(df35$y), rep(NA,7)),state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.family ="student", #A Bayesian Analysis of Time-Series Event Count Dataniter =3e4, #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.seed=2125)## Since we have predictor variables in the model, we need to explicitly# make their coefficients time-varying using AddDynamicRegression().impact3d_ratio_resp <-CausalImpact(bsts.model = model2d_ratio,model.args =list(prior.level.sd=.1, dynamic.regression=T),post.period.response = df35_post$y)plot(impact3d_ratio_resp, "original")
# Model 2ss2a <-list()# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TRENDss2a <-AddLocalLevel(ss2a,c(as.numeric(df33$y), rep(NA,7)) ) ## Add weekly seasonalss2a <-AddSeasonal(ss2a,c(as.numeric(df33$y), rep(NA,7)), nseasons=7, season.duration =1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑOss2a <-AddAutoAr(ss2a, y =c(as.numeric(df33$y), rep(NA,7)), lags =10) #NO PUEDO AREGAR AR1 CON POISSON# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).# Add external regressors explicitly (e.g., temperatures, precipitation)covariates33 <-cbind(# y = as.numeric(df35$y),mean_min_temp_mean_lin = df33$mean_min_temp_mean_lin,mean_max_temp_mean_lin = df33$mean_max_temp_mean_lin,mean_median_lag_2_prec_median_imp = df33$mean_median_lag_2_prec_median_imp,mean_median_total_precip_median_lin= df33$mean_median_total_precip_median_lin,Y_lag10 = df33$Y_lag10)# Ensure names explicitlycolnames(covariates33) <-c(#"y","mean_min_temp_mean_lin","mean_max_temp_mean_lin","mean_median_lag_2_prec_median_imp","mean_median_total_precip_median_lin","Y_lag10")covariates33<-as.data.frame(covariates33)covariates_future33 <-as.matrix(df33_post[1:7, c("mean_min_temp_mean_lin","mean_max_temp_mean_lin","mean_median_lag_2_prec_median_imp","mean_median_total_precip_median_lin","Y_lag10")])covariates_full33 <-rbind(covariates33, covariates_future33)# Include regression termsss2a <-AddDynamicRegression(ss2a, covariates_full33)model2a_ratio <-bsts(c(as.numeric(df33$y), rep(NA,7)),state.specification = ss2a, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.family ="student", #A Bayesian Analysis of Time-Series Event Count Dataniter =3e4, #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.seed=2125)## Since we have predictor variables in the model, we need to explicitly# make their coefficients time-varying using AddDynamicRegression().impact3a_ratio_resp <-CausalImpact(bsts.model = model2a_ratio,model.args =list(prior.level.sd=.1, dynamic.regression=T),post.period.response = df33_post$y[1:7])plot(impact3a_ratio_resp, "original")
# Model 2ss2ab <-list()# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TRENDss2a <-AddLocalLevel(ss2ab,c(as.numeric(df33$y), rep(NA,3)) ) ## Add weekly seasonalss2ab <-AddSeasonal(ss2ab,c(as.numeric(df33$y), rep(NA,3)), nseasons=7, season.duration =1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑOss2ab <-AddAutoAr(ss2ab, y =c(as.numeric(df33$y), rep(NA,3)), lags =10) #NO PUEDO AREGAR AR1 CON POISSON# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).covariates_future33b <-as.matrix(df33_post[1:3, c("mean_min_temp_mean_lin","mean_max_temp_mean_lin","mean_median_lag_2_prec_median_imp","mean_median_total_precip_median_lin","Y_lag10")])covariates_full33b <-rbind(covariates33, covariates_future33b)# Include regression termsss2ab <-AddDynamicRegression(ss2ab, covariates_full33b)model2ab_ratio <-bsts(c(as.numeric(df33$y), rep(NA,3)),state.specification = ss2ab, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.family ="student", #A Bayesian Analysis of Time-Series Event Count Dataniter =3e4, #burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.seed=2125)## Since we have predictor variables in the model, we need to explicitly# make their coefficients time-varying using AddDynamicRegression().impact3ab_ratio_resp <-CausalImpact(bsts.model = model2ab_ratio,model.args =list(prior.level.sd=.1, dynamic.regression=T),post.period.response = df33_post$y[1:3])plot(impact3ab_ratio_resp, "original")
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact33b_qc' no encontrado
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': objeto 'impact2_33b_qc' no encontrado
Brodersen KH, Gallusser F, Koehler J, Remy N, Scott SL. Inferring causal impact using bayesian structural time-series models. The Annals of Applied Statistics. 2015;9(1). doi:10.1214/14-aoas788.